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Dynamical facilitation theory assumes short-ranged dynamical constraints to be the essential 
feature of supercooled liquids and draws much of its conclusions from the study of kinetically con- 
strained models. While deceptively simple, these models predict the existence of trajectories that 
maintain a high overlap with their initial state over many structural relaxation times. We use 
molecular dynamics simulations combined with importance sampling in trajectory space to test 
this prediction through counting long-lived particle displacements. For observation times longer 
than the structural relaxation time exponential tails emerge in the probability distribution of this 
number. Reweighting trajectories towards low mobility corresponds to a phase transition into an 
inactive phase. While dynamics in these two phases is drastically different structural measures show 
only slight differences. We discuss the choice of dynamical order parameter and give a possible 
explanation for the microscopic origin of the effective dynamical constraints. 



Since crystallization occurs through nucleation virtu- 
ally any liquid can be supercooled below its melting tem- 
perature. But some liquids never become crystals. Their 
viscosity increases dramatically and at some point inter- 
nal relaxation cannot keep up with the cooling anymore 
and they fall out of equilibrium, reaching a state we call 
a glass (for reviews see Refs. [HH]). While in principle 
protocol-dependent, the temperature range at which the 
transition occurs is narrow and nearly a material prop- 
erty. The idea that this glass transition is a, or deter- 
mined by a, thermodynamic transition has influenced the 
theoretical studies of glasses for decades. However, ex- 
perimentally determined structure factors of supercooled 
liquids show little to no change while approaching the 
glass transition. If not global structure, what is the ori- 
gin of slow dynamics? 

The arguably most striking feature of supercooled liq- 
uids is the emergence of dynamical heterogeneity (see 
Ref. [3] for a review) below an onset temperature, i.e., 
while large regions of the liquid are jammed structural 
relaxation continues through regions which are less rigid. 
This phenomenon has been observed directly in, e.g., col- 
loidal glasses [1] and granular systems [S]. While simple 
liquids above the onset temperature are well described by 
a mean- field theory [6j , dynamical heterogeneity leads to 
different environments for different particles. This ob- 
servation forms the foundation for dynamical facilitation 
theory [7] , a theory of the glass transition as a dynami- 
cal phenomenon, see Ref. [8] for a review and references. 
While the interplay between structure and particle dy- 
namics is complicated, the glass transition is indepen- 
dent of much of these details and dynamics is dominated 
by effective constraints restricting the accessible regions 
in space-time. The glass transition is controled by the 
number of excitations marking locally weak or soft re- 
gions able to reorganize, and not by a thermodynamic 
variable. 

Crucial for dynamical facilitation theory is the notion 
of a mobility field coarse-grained both in time and space. 



A convenient caricature of mobility fields are spin-like 
excitations on a lattice [9 . The effect of the crowded 
environment on particle motion is mimicked by kinetic 
constraints, i.e., for a spin to change its state it must be 
facilitated by one (or more) neighboring excited spin(s). 
These dynamical rules suffice to give rise to dynamical 
heterogeneity and a dramatic slow-down of relaxation in 
the absence of thermodynamic transitions. 

Structural relaxation is the decorrelation of particle 
positions with their initial positions over time. In a 
dense liquid particle motion is strongly hindered by sur- 
rounding particles leading mainly to vibrational motion, 
short-lived excursions, and rare, collective, long-lived 
particle displacements that could be described as "cage 
breaks" [TO]. One approach to glassy dynamics is to find 
strategies to predict long-time motion from short-time vi- 
brations probing the local structure [HI HI]- In this paper 
we pursue a different route and focus on the long-lived 
displacements as recorders of excitations in the sense of 
the simple lattice models. We introduce a binary mobil- 
ity field but use molecular dynamics to determine its time 
evolution instead of postulated, and necessarily idealized, 
dynamics. 

The paper is organized as follows: In Sec. [T] we give 
a brief reminder on kinetically constrained models. In 
Sec.|ll]we combine the technique introduced in Ref. |13j 
to record excitations in atomistic models of glass form- 
ers with the fundamental ideas outlined in Refs. [13] and 
|15j . We thus confirm predictions made previously by the 
study of kinetically constrained lattice models. In partic- 
ular, for observation times much longer than the relax- 
ation time we find exponential tails in the distribution 
of mobile particles and a phase transition upon varying 
a field that couples to mobility. Such dynamical phase 
transitions [16] have been studied analytically and nu- 
merically for kinetic constrained lattice models [T7| (TB] , 
spin-glasses [19], and have also been found in atomistic 
glass formers [20]. In Sec. Ill we then study some struc- 



tural and dynamical properties of the inactive phase in 
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more detail before coming to the conclusions in Sec. IV 



I. KINETICALLY CONSTRAINED LATTICE 
MODELS 

We start with a brief reminder on two popular vari- 
ants of kinetically constrained lattice models [H]: the 
Frederickson- Andersen (FA) |9j and the East model [22] . 
Both models have a trivial energy function E({ni}) = 
J^jTi; where J sets the energy scale and ni is either 
1 or indicating whether site I is excited or not. Ex- 
citations correspond to regions of the supercooled liquid 
where particles are unjammed and mobile. The effects 
of local jamming are incorporated by kinetic constraints: 
in the FA model a site can change state only if a neigh- 
boring site is excited, i.e., ni±i ^ 1; in the East model 
the site / -I- 1 must be excited. Evolution of the mobility 
field {ni} is assumed to be Markovian. A simple choice 
for the rates is the following. Annihilation ri; : 1 of 
an excitation occurs with rate 1 defining the time scale. 
Detailed balance then implies the rate e""^/"^ for the cre- 
ation of an excitation (if the constraint is fulfilled) . The 
equilibrium concentration of excitations is 
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with respect to the reduced temperature 1/T = 1 /T — 
1/Tq. The temperature Tq marks the onset of dynamical 
heterogeneity. 

As a result of constraints, excitations move in lines. 
Fluctuations in excitations lead to coalescing and branch- 
ing of excitation lines. These changes of state, or "kinks" , 
provide the means to detect excitations through particle 
motion as discussed in the next section. However, the ab- 
sence of detectable motion does not signify the absence 
of excitations which, without other excitations reaching 
out, are quiescent as a direct consequence of the kinetic 
constraints. The phase transition that we detail in this 
paper is a transition between two phases of markedly dif- 
ferent concentrations of kinks or fluctuations. 



II. PARTICLE MOBILITY AS SPACE-TIME 
ORDER PARAMETER 

A. Excited particles 

We have performed extensive molecular dynamics sim- 
ulations on the Kob- Anderson binary mixture ^23j, a 
popular model for atomistic glass formers [211 di] (see 
appendix IA] for details). It is composed of 80% large 
(A) and 20% small (B) particles. Simulations are run at 
temperature T = 0.6 well below the onset temperature 
To ~ 0.88 of heterogeneous dynamics. For comparison, 
T ~ 0.435 is an estimate of the glass transition tempera- 
ture for this system [521 US] . We have chosen the higher 
temperature and therefore only moderately supercooled 



state point to be able to run millions of trajectories over 
a couple of structural relaxation times. For that state 
point, the structural relaxation time is ~ 24.5 as mea- 
sured by the decay F{Ta) = 1/e of the intermediate scat- 
tering function (see Fig. [I^) 
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Particle motion is measured on the length scale 27r/|q| ~ 
1.058 corresponding to the peak position of the pair dis- 
tribution function. Fig. [T^ demonstrates that at T = 0.6 
the structural relaxation for the system sizes considered 
here shows no finite size effects (see also Ref. [26]). 

Kinetically constrained models are minimal models in- 
corporating the crucial ingredient of hindered motion in 
dense, nearly jammed liquids. What they did not provide 
so far is a concrete prescription on how to map a set of 
particle positions {r;} onto a binary mobility field {n;}. 
To establish such a mapping we follow Ref. J3j and use 
long-lived particle displacements of a given length a as 
recorders of excitations. To this end we define the single 
particle indicator function 



hi{t) = Q{\f i{t) - f lit -At)\~ a), 
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where 0{x) is the unit step function and /i; = 1 if the 
Ith particle has moved further than the distance a in the 
time interval At and hi — otherwise. In the following 
we will call particles with hi = 1 mobile, or excited to 
emphasize their role as recorders of underlying excita- 
tions. To distinguish non-trivial particle displacements 
from mere vibrations we use the inherent structure posi- 
tions {f;} [27]. The inherent structure of a configuration 
is obtained by steepest descent to the nearest minimum 
in the potential energy landscape. Even though there is 
a hierarchy of motion on different length scales [13] , here 
we focus on the single length a — 0.3. The commitment 
time At = 1.5 is then chosen to be large enough so that a 
particle can commit to a new position but small enough 
so that we do not count multiple jumps. 

The instantaneous density of mobile particles is 
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with mean c = (c) = (hi). In Fig.[T|D we demonstrate that 
the temperature dependence of this mean density indeed 
follows the prediction for excitations Eq. ([!]). The fitted 
onset temperature Tq ~ 0.88 agrees excellently with pre- 
vious estimates [T3l|28]. The fitted energy is J ~ 3.9 j49| . 



B. Low mobility tails 

Our objects of interest are trajectories X = {r;(t) : 
^ t ^ tobs} of fixed length tobs- To quantify the 
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FIG. 1: (a) Intermediate scattering function Eq. ([2| at T = 0.6 for two system sizes displaying the two-step decay typical for 
supercooled liquids, (b) The mean density of excited particles c (symbols) as a function of inverse temperature 1/T below the 
onset temperature To ~ 0.88. The dashed line is the fit to Eq. M. (c+d) Logarithm of the probability p(c) of the intensive 
order parameter c scaled by the particle number for (c) A'^ = 216 and (d) = 343 particles, and different observation times 
tohs- The dashed lines with slope 7 = 0.58 indicate the exponential tails. 



amount of mobility in a trajectory we count excited par- 
ticles through the order parameter 
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with equally spaced ti = ilS.t and observation time 
^obs = Kl\t. In Fig. [1] we show distributions p{c) for 
different system sizes and different observation times ob- 
tained through umbrella sampling combined with replica 
exchange (see appendix [b|). For trajectories in which mo- 
tion decorrelates on a time scale much shorter than tobs 
we would sum over many independent events. Follow- 
ing the central limit theorem the probability distribution 
p(c) then approaches a Gaussian. Indeed, for moder- 
ate observation times Tq < tobs ^ ^obs l^^rger than the 
relaxation time Tq, but below a cross-over time t*^^^ we 
find such Gaussian distributions (see iobs = 75 w Sr^ in 
Fig. [1]). Increasing tobs ^ ^obs observe two effects: 
for small c exponential tails emerge and the shape of 
p{c) becomes non-concave. The physical picture is that 
highly constrained dynamics facilitates the creation of ex- 
tended immobile regions, i.e., compared to uncorrelated 
dynamics it is easier to "remove" mobility from a given 
space-time volume. 

The explanation why the tails of p{c) are exponential 
is as follows P^. Assuming that excitations are non- 
interacting the probability to find at t = an immobile 
region of i particles is proportional to e~^'^^ with a geo- 
metric factor 7 independent of N and tobs- Combining 
this probability with C ~ K{N — i)cior the case that this 
"bubble" persists leads to Inp(c) — jNc plus an offset. 
In Ref. |13j , evidence is presented that the assumption of 
non-interacting excitations is indeed a good approxima- 
tion. Of course, not all bubbles span the entire trajectory 
connecting the initial with the final state. The temporal 
extent t*,^g of a typical bubble grows proportionally to 
the mean persistence time, i.e., the mean time a particle 
remains at its initial inherent structure position. 



The order parameter Eq. ([s]) is purely dynamical. 
Since in a crystal particle mobility is low and motion 
restricted to defects such an order parameter cannot dis- 
criminate between a low activity amorphous and a low 
activity crystalline phase. Therefore, we also monitor 
the structure through the orientational order parameter 
■06(0 defined in Eq. (Dl). To allow local fluctuations in 



structure but prevent global long-range order we use the 
value of 
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averaged over the Na — Q.8N large particles of the final 
configuration of trajectories. We reject all trajectories 
with > 0.45. (see also Fig.jsJ; below) 



C. Active— inactive phase transition 

So far we have established that trajectories contribut- 
ing to the exponential tails in Fig. [Tjare those which re- 
member their initial conditions and do not relax within 
the observation time tobs- At least formally and in com- 
puter simulations we can apply a bias in the space of tra- 
jectories to stabilize these low mobility trajectories. The 
fact that the distributions p{c) are non-concave implies a 
phase transition between an active phase corresponding 
to the liquid in which motion is plentiful, and an inactive 
phase of low particle mobility. To provide a link with 
traditional thermodynamics, and the Ising model in par- 
ticular, imagine the {hi} to be spins on a lattice with 
the order parameter C taking the role of the magneti- 
zation [291 . Below the critical temperature the system 
undergoes a first-order transition between a disordered 
phase of low magnetization and an ordered phase of high 
magnetization through applying a field (which we will 
call s in the following). While the statistical treatment 
is analogous, the underlying physics of a supercooled liq- 
uid is of course different from the Ising model. In our 
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FIG. 2: Mean fraction of excited particles (top) and suscepti- 
bility (bottom) vs. the biasing field s for selected observation 
times (from right to left: fobs ~ 300, 375, 600) and system 
sizes N = 216 (solid lines) and A'^ — 343 (dashed lines). For 
clarity we only show error bars for the peak values of the sus- 
ceptibility x{^)- Upper inset: Dependence of the coexistence 
field s* on trajectory length iobs and system size A'^. The 
dashed lines are fits to s* = 8% + a/tobs- The fit parameters 
are S216 — 5.1 x 10~^ and S343 ~ 5.8 x 10~^. From these 
results we estimate the coexistence field to be sj^ ~ 5 x lO"** 
(arrow) in the limit fobs — >■ 00. Lower inset: System size 
dependence of x* ■ 



case the lattice extends over space and time, and the in- 
teractions between "spins" is due to short-ranged forces, 
geometrical confinement, and the thereof resulting tem- 
poral correlations of particle motion. 

In Fig. [2] we plot the mean fraction of excited particles 
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and the susceptibility x(s) = —d{c)s/ds vs. the bias- 
ing field s. The plot shows that the density of mo- 
bile particles drops from c ~ 0.11 at s = to about 
Cin — 0.01 for s ^ s* . For small s we can expand the 
mean (c)^ ^ c- ns with n = [(C^) - {Cf]/{NK). The 
linear behavior around s = in Fig. [2] therefore, reflects 
the Gaussian nature of the liquid phase. The coexistence 
field s* is obtained from the peak position of x(s) maxi- 
mizing the fluctuations of the order parameter. Increas- 
ing either the number of particles N or the observation 
time fobs sharpens the transition. However, space and 
time are not symmetric. At least in part, this asymme- 
try reflects that we employ periodic boundary conditions 
in space whereas trajectories can have quite distinct ini- 
tial and final states. This leads to temporal boundary 
effects enhancing the mobility at the beginning and the 
end of the trajectory [16]. For fixed TV one expects to 
leading order s* = s*pf + ©(l/tobs) as shown in the upper 
inset of Fig. [2] From the fits we estimate the limiting 



coexistence field « 5 x 10~^ to be small but nonzero. 
Finally, in the lower inset of Fig. [2] we demonstrate the 
finite-size scaling of the peak values x* = x(s*) plotted 
vs. the space-time volume Ntohs — NKAt. The dashed 
line corresponds to x* ~ tobs(0.065 + 1.08 x lO^^A^fobs), 
which suggests a first-order transition with a diverging 
susceptibility and a discontinuous jump of (c)^ at in 
the limit of large N and/or fobs- 



D. Choice of order parameter 

To check whether the transition into an inactive phase 
is robust with respect to the way we measure mobility in 
trajectories, we have also considered the dynamical order 
parameter 

K N 
i=l 1 = 1 

(8) 

This order parameter was used in Ref. |20| to demon- 
strate for the first time a transition between a high ac- 
tivity and a low activity phase in an atomistic model. 
It sums over the short-time mean-square displacements 
of particles, which obscures the separation of vibrations 
from reorganization events that lead to structural relax- 
ation. Nevertheless, as demonstrated in Fig.jSj an abrupt 
transition from high to low activity can still be observed. 
To be specific, we define a new ensemble 



-cr/C\ 
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where we denote the biasing field coupling to /C by cr 
and A is any observable. As an illustration for N — 216 
and fobs = 450 the mean {k)cr is plotted in Fig. [sj It 
resembles the curves shown in Fig. [2] and drops abruptly 
around cr* ~ 0.0077. Moreover, we find exponential tails 




FIG. 3: Mean intensive activity {fc)o- vs. the biasing field a 
for A'' — 216 and fobs = 450 (solid line). For comparision, 
the mean fraction of excited particles (right axis) is shown for 
both the ensembles defined through IC {{c)cr, dashed line) and 
C {{c)s, dash-dotted line). Inset: Probability distribution of 
k. The dashed line indicates the exponential tail. 
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in the probability distribution of k plotted in the inset 
of Fig. [sj In addition to (fc)o- we determine the density 
of mobile particles (c)o-, which closely follows the for- 
mer. However, while the activity drops by a factor of 
less than two, the number of excited particles is reduced 
by more than one order of magnitude. The comparision 
with the curve (c)s obtained in the C-ensemble shows 
that the mean fraction of mobile particles approaches the 
same value (within uncertainties) in both ensembles. We, 
therefore, conclude that both measures C and /C prepare 
the same inactive phase through applying a biasing field 
in trajectory space. The transition is more pronounced 
in C since vibrational motion, which does not cease in the 
inactive phase evolving at the same temperature as the 
active phase, contributes significantly to Eq. (Isl). 
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III. PROPERTIES OF THE INACTIVE PHASE 
A. Nucleation of activity 

An important consequence of first order phase transi- 
tions is nucleation: a system crossing the transition line, 
although the new phase is favored, has to pay a penalty 
for interfaces and one has to wait for a large enough nu- 
cleus to appear spontaneously. Translated to the present 
case one might ask what happens if at time t = tobs for 
s ^ s* we turn off the field s. In the picture of facili- 
tated dynamics a new excitation can only appear close 
to an existing excitation. Since both density of excita- 
tions and kinks, as recorded through (c)s, are drastically 
reduced in the inactive phase we have to wait a certain 
time before excitations "percolate" through the system 
and it returns to the liquid phase. The nucleation of ac- 
tivity, therefore, is conceptually different from, say, the 
nucleation of a crystal which is determined by a single 
large barrier. 

To demonstrate this "melting" of jammed configura- 
tions we prepare trajectories at s = 0.01. From these 
trajectories we then take a single configuration, random- 
ize velocities, and run 10,000 unbiased trajectories out of 
this initial configuration (see also the isoconfigurational 
ensemble [30] )■ In Fig. |4ji a single trajectory is shown. 
Clearly, the system remains inactive for many structural 
relaxation times (in this example — 15tq,) as measured 
by c and then, suddenly, becomes active again with large 
fluctuations of c. In Fig. |4]3 a different trajectory out of 
the same initial configuration shows a more gradual tran- 
sition from inactive to active. In Fig. [4]: and d we show 
distributions of waiting times for two different initial 
states. These distributions are clearly non-exponential, 
which is consistent with a variable step process as exci- 
tations reach out and reconnect. A detailed comparision 
of these distributions to predictions from kinetically con- 
strained models is left for a future study. 



FIG. 4: Melting of the inactive phase: (a) Fraction of mobile 
particles c{t) vs. time for an initial configuration prepared 
with the s-ensemble at s = 0.01. The systems remains in- 
active up to time tw (arrow), when it suddenly unjams. The 
dashed line indicates the average density c. (b) Another melt- 
ing trajectory showing a more gradual "thawing", (c) The 
distribution of tw for 10,000 trajectories starting out of the 
same initial configuration, (d) A different initial configura- 
tion showing a much broader distribution of waiting times. 



B. Local structure 

In the introduction we have emphasized that global 
structural differences between liquid and glass are at 
most minuscule. However, the fact that the system re- 
mains jammed for particle configurations taken from tra- 
jectories prepared at s ^ s* indicates that there is a 
structural difference between these configurations and 
configurations typically visited in the liquid phase. To 
make this more quantitative we sample trajectories at 
fixed s and compare the structures as measured by three 
different methods, see Fig. [5] The pair distribution func- 
tion gAAir) for the large (A) particles shown in Fig. [5^ 
demonstrates that liquid and inactive phase are globally 
indistinguishable. Small differences are seen for the small 
(B) particles in Fig. ^jp. Beyond the simple two-point 
functions we also consider the histogram of the bond- 
order parameter tpQ as defined in appendix |D] plotted in 
Fig. [5];. This order parameter is a convenient measure 
for long-range order. For every particle it quantifies its 
local order with ip^ — 1 for a particle in a perfect crys- 
tal. All measures clearly show that the inactive phase is 
amorphous. 

In Fig. [5ji we show that the potential energy per par- 
ticle and the density of mobile particles are uncorrelated 
in both phases. While in the inactive phase the potential 
energy of particles is typically lower, the mean differ- 
ence « O.l is much less than the vibrational contribution 
~ 1.5T separating real space potential energies from the 
inherent state energies. Moreover, as demonstrated in 
Fig. [5}l, there is still an overlap of potential energies be- 
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FIG. 5: Different measures of the structure in the hquid (s = 0) and the inactive phase (s = 0.01) for A'' = 216 and 
iobs = 600: (a) Radial distribution function for the large A particles and (b) for the small B particles. The arrow marks the 
peak corresponding to a B-B bond with 3 common A neighbors, see sketch and main text, (c) Distribution of the structural 
order parameter i/jg measuring long range order. A particle in a perfect crystal would have V'e = 1- (d) Scatter plot of potential 
energy (PE) per particle versus the concentration of excited particles c for both actual positions and inherent states. Red points 
are from the ensemble of active states and black points are from the ensemble of inactive states. The harmonic contribution 
3r/2 to the potential energy is indicated. 



tween both phases. Hence, we conclude that particles are 
not trapped energetically but rather due to geometrical 
constraints. 

Differences in structure are picked up by the pair dis- 
tribution function for the small (B) particles plotted in 
Fig. [sja. It has been shown that the peaks of gBB{r) for 
binary mixtures can be assigned to certain local struc- 
tures: two bonded B particles sharing m common A 
neighbors SJ. • Of particular interest is the second peak 
corresponding to m = 3 since it indicates icosahedral 
coordination shells. Fig. [SJd shows that in the inactive 
phase this local structure occurs more often compared 
to s = 0. This is consistent with recent observations 
of short-ranged structures in supercooled binary mix- 
tures [321 [33] . Slow relaxation is attributed to reorgani- 
zation of particles bound in these structures. Moreover, 
the drop of ~ 0.1 in the potential energy of inherent 
states agrees quantitatively with the drop associated to 
the formation of these structures [33] (albeit for a slightly 
different model). 



is wi = 1. We follow a single particle along a trajectory 
and through 
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(11) 

we count the number of excited particles that have been 
created in a sphere with radius r around the tagged par- 
ticle under the condition that the tagged particle itself 
had been excited in the preceding time slice. We define 
the transfer function 
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The ratio ^(s)//i(0) is plotted in the inset of Fig. [6] us- 
ing r = 1.5 roughly corresponding to the first coordina- 
tion shell. It shows that the probability that a particle 
becomes excited close to an already mobile particle in- 
creases in the inactive phase. 



C. Dynamical facilitation 

We finally study the behavior of facilitation when go- 
ing from the active to the inactive phase. First, we note 
that the fraction of excited particles as plotted in Fig.[6|is 
independent of temperature in the inactive phase. This 
indicates that the dynamics in the inactive phase is de- 
coupled from the externally fixed temperature. Second, 
we study the degree to which particle motion is facili- 
tated. Different methods have been reported in the lit- 
erature including a mobility transfer function |34| and 
the facilitation volume [13j . In the spirit of a mobility 
transfer we consider the set of newly excited particles for 
which the binary indicator function 
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FIG. 6: Mean fraction of excited particles (c)s at three differ- 
ent temperatures T. Inset: The normalized transfer function 
Eq. (121 showing an increase of facilitation at large s. (all 
data for A^ = 216 and ioba = 300) 
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Putting all observations together the following picture 
of dynamics in the inactive phase emerges: The persis- 
tence time exceeds the observation time and most par- 
ticles maintain a high overlap with their initial position. 
However, some activity continues in isolated regions. The 
fraction of these mobile particles is decoupled from the 
temperature. Particles do not become mobile (or im- 
mobile) at random but are facilitated through existing 
mobile particles in their vicinity. Turning off the s-field 
these remaining mobile particles are the seeds from which 
excitations can reconnect before the system returns to its 
fluid state. 



IV. CONCLUSIONS 

The separation between fast inter-basin vibrations and 
slow, activated transitions between inherent structures 
(or meta-basins [3S]) is the essence of the energy land- 
scape paradigm [5J . It implies a time evolution that 
is dominated by rare thermal fluctuations that carry the 
system from one minimum over a barrier into a neighbor- 
ing minium |37) . However, there is mounting evidence 
that such a mechanism competes with, or is even shad- 
owed by, relaxation that occurs through channels that 
present only low energetic barriers but which are rare 
and found through "surging" particle motion: examples 
are string-like motion ^38j and participation maps of low- 
frequency modes [30l |39] . Dynamical facilitation theory 
postulates "excitations" to be the fundamental objects 
describing such structural weaknesses, and facilitation to 
be the dominant mechanism at low concentrations of ex- 
citations. 

The energy landscape picture assumes that the way the 
landscape is sampled is governed by temperature. From 
the evidence presented here, we arrive at a seemingly 
different perspective: the concentration of excitations de- 
termines relaxation. At constant temperature the system 
can be forced into an inactive glassy phase through either 
removing excitations or suppressing fluctuations leading 
to "frozen" excitations that are quiescent. In this inac- 
tive phase particles vibrate around local energy minima, 
the statistics of which is consistent with the externally 
fixed temperature. In contrast, transitions between lo- 
cal minima, or inherent states, are rare and decoupled 
from this temperature. It appears that these jammed 
states, in which excitations are arrested, can be created 
rather easily through local particle rearrangements that 
do not affect the global structure. In Ref. it is shown 
that these states are mechanically more stable than fluid 
states at a lower temperature. Here we have demon- 
strated that the melting of jammed states is not consis- 
tent with a single crossing of a large free energy barrier 
but that it is rather a multi-step process. We attribute 
this multi-step process to the "unfreezing" of excitations, 
an interpretation that is supported by an increased de- 
gree of facilitation in these jammed states. These obser- 
vations, together with the emergence of exponential tails 



equivalent to those observed in kinetically constrained 
models, leads us to the conclusion that the transition 
is indeed caused by local dynamic constraints. The pre- 
cise pathways and microscopic mechanisms of the particle 
rearrangements underlying the active-inactive transition 
are left for future studies. 

As a final note we emphasize that the active to inactive 
transition is reminiscent of the glass transition. The fun- 
damental difference is that the transition demonstrated 
in this paper is a transition controlled by a field coupling 
to a dynamical observable, while the experimental glass 
transition is controlled by rate of temperature decrease. 
The connection between the two remains to be quanti- 
fied. 
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Appendix A: Simulation details 

We have performed extensive molecular dynamics sim- 
ulations on the binary mixture p3]. It is composed of 
80% large (A) and 20% small (B) particles possessing the 
same mass m. Particles interact through the continuous 
truncated and shifted Lennard- Jones potentials 

{"Lj(f';ea/3,0-a/3)- 
MLj(2.5cra^;ea^,Cra/3), T < 2.5crQ,3 
0, r> 2.b<Jaf} 

where MLj(r;£,CT) — Ae[[a/rY'^ — (cr/r)^]. The parame- 
ters read ctaa = cr, ctab = O.Scr, cbb — O.SScr, £aa — £, 
Sab — 1.5e, and £bb — O.Se. Simulations are run at 
constant volume V , constant temperature T, and con- 
stant number density N/V = 1.2 with N the number 
of particles at positions {r;}. Throughout the paper we 
employ reduced Lennard- Jones units with respect to the 
large particles, i.e., we measure length in units of cr, en- 
ergy in units of e, time in units of -^/mcr^/e, and we set 
Boltzmann's constant to unity. 

Newton's equations of motion are integrated through 
the velocity Verlet algorithm with time step 0.005 using 
LAMMPS |41j . For energy minimization we employ the 
FIRE algorithm gg. 



Appendix B: Importance sampling 

Just as standard Monte Carlo simulations 43J do im- 
portance sampling of configurations, the methods we em- 
ploy do importance sampling of trajectories. Specifically, 
we harvest new trajectories using moves from transition 
path sampling [44l [45] . These moves preserve the equilib- 
rium weight -Po[^] of trajectories. Hence, accepting or re- 
jecting a move X ^ X' according to the usual Metropolis 
criterion 



min < 1 



^-w{C[X']) 

e-^ic[x]) 



generates an ensemble of trajectories with weight 

P[X] - Po[X]e'""^'^™ . 

Here, C[X] is a dynamical order parameter that calculates 
a real number from trajectory X, see Eqs. ([5| and (|8]); 
and 'w{x) is the weight function. 

For a single trajectory we store K + I configurations 
at times U = iAt with i = 0, . . . ,K. We employ the 
"massive stochastic collision" thermostat, i.e., all veloci- 
ties are randomized at these times and the center-of-mass 
velocity is subtracted. The use of a stochastic thermo- 
stat allows us to perform transition path sampling with 
so-called 'half moves' as described in detail in Ref. 
In order to efficiently sample probability distributions of 
the order parameter we use the quadratic form 

w{x)='^{x-xof. (Bl) 

To speed up the sampling of trajectory space and de- 
crease the correlations of subsequently generated trajec- 
tories we use replica exchange between iVrcp = 8 replicas 
with different values of {xq}. Hence, a single cycle con- 
sists of two consecutive steps: (i) every replica generates 
a new trajectory which is either accepted (and replaces 
the previous trajectory) or rejected, and (ii) trajectories 
are swapped between replicas. To obtain good mixing 
we attempt 8^ swaps between all of the replicas, not only 
neighbors. Data has been acquired from two indepen- 
dent runs, i.e., two independent seed trajectories (except 
K = 50, which is from one run). For a single run we 
let the trajectories relax for cycles and then recorded 
A^t trajectories. Table [l] shows an overview for the data 
gathered to produce Figs, [l] and [2j 

Appendix C: MBAR 

To calculate distributions and expectation values from 
raw data we use Shirts' and Chodera's multistate Bennett 
acceptance ratio (MBAR) method [15] and its extension 
to path ensembles [47j . For completeness we briefly sum- 
marize the method. We solve 



TABLE I: System sizes studied and number of harvested tra- 
jectories. 
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"Low activity umbrellas might not have fully equilibrated. 

self-consistently for the set of "free energies" {fi}- Here, 
Wi{x) is the weight function corresponding to replica i, 
and Xjn = Wj{C[Xn\) is the value of the order parameter 
along the nth trajectory of replica j. While in principle 
MBAR provides an error estimate it requires indepen- 
dent samples, whereas consecutive trajectories obtained 
using the method described in the previous section are 
highly correlated. The statistical inefficiency of replica j 



is gj = 1 + 2Tj , where Tj is the correlation time (in unit 
samples) of some representative observable (here we use 
c). One possibility to obtain independent samples is to 
subsample the data series with stride gj . Here we use all 
data but weigh replicas according to their relative sta- 
tistical inefficiencies. Errors are estimated by splitting 
the data into chunks of A't = 10, 000 trajectories and 
calculating the standard error of expectation values. 

Expectation values of an observable A are calculated 
through 



,-[w{x,„)-f] 



Here, w(a;) and / can correspond to one of the repli- 
cas, i.e., w — Wi and f = fi. The advantage of MBAR 
is that we can employ an in principle arbitrary weight 
function w{x) (given sufhcient statistical weight e"™*^^-* 
of the sampled data) with 



N„ 



Nt 



(CI) 



Probabilities pi = (xi) are calculated as the expectation 
value of an indicator function Xi{^) which is 1 if the value 
X falls into bin i and otherwise. In particular, the 
distributions shown in Fig. [T] correspond to the unbiased 
ensemble -Po[^] with w{x) — 0. 

The curves shown in Fig. [2] for the mean value (c)s are 
obtained through using the weight function w(x) = sx in 
Eq. (CI). In order to sample trajectories at fixed s we 



use this weight function instead of Eq. (Bl) for a chain 
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of replicas with different s values ranging from s = to 
s = 0.01. To obtain a set of independent trajectories we 
keep only every lOOOtli trajectory for analysis. 

Appendix D: Orientational order 

To quantify orientational order we follow Ref. [IS]. For 
each particle k a complex vector 




is defined, where Yim are spherical harmonics and the 
angles 9kk' and <j)kk' describe the orientation of the dis- 
placement vector between particles k and k' with respect 



to a fixed reference frame. The sum is over all iVb neigh- 
bors in the first coordination shell with radius 1.42. The 
normalized scalar product of the g- vectors is 

The average over neighbors 

^ k' = l 

is the bond order parameter. It is '06 = 1 for a particle in 
a perfect crystal and acquires a broad distribution with 
mean 0.2 — 0.3 for particles in a disordered environment. 
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